####################
####################
###Moran's I Test###
####################
####################
library(doBy)
library(spdep)
library(sf)
library(dplyr)
library(lmtest)
library(sandwich)

# Define your subsets
subsets <- list(
  ag = subset(med.dat, ag_crops > 0),
  forest = subset(med.dat, Forest.coverage_bin > 0),
  other = subset(med.dat, Forest.coverage_bin == 0 | ag_crops == 0)
)

# Prepare to store results
results <- data.frame(subset_name = character(),
                      morans_I = numeric(),
                      morans_I_se = numeric(),
                      mean_confirm = numeric(),
                      p_value = numeric(),
                      significant = character(),
                      stringsAsFactors = FALSE)

# Load shapefile once
polygons.prio <- st_read("priogrid_cellshp/priogrid_cell.shp", quiet = TRUE) %>%
  dplyr::select(gid)

# Loop through subsets
for (name in names(subsets)) {
  cat("Processing subset:", name, "\n")
  
  subset_df <- subsets[[name]]
  
  # Collapse to one value per gid (e.g., total outbreaks)
  ag_df <- summaryBy(clear_zoon_confirm ~ gid, data = subset_df, FUN = sum, keep.names = TRUE, na.rm = TRUE)
  
  # Subset the shapefile to the GIDs in this subset
  polygons.slice <- subset(polygons.prio, gid %in% ag_df$gid)
  rownames(polygons.slice) <- 1:nrow(polygons.slice)
  
  if (nrow(polygons.slice) == 0) {
    cat("No polygons matched for subset", name, "- skipping.\n")
    next
  }
  
  # Create spatial weights
  slice.sp <- as(polygons.slice, "Spatial")
  slice.nb <- poly2nb(slice.sp, queen = TRUE)
  slice.listw <- nb2listw(slice.nb, style = "W", zero.policy = TRUE)
  
  # Match order of ag_df to polygon order (just in case)
  ag_df <- ag_df[match(polygons.slice$gid, ag_df$gid), ]
  
  # Spatial lag
  inc.lag <- lag.listw(slice.listw, ag_df$clear_zoon_confirm)
  
  # Plot
  jpeg(paste0("morans_", name, ".jpeg"), width = 6, height = 6, units = 'in', res = 500)
  plot(inc.lag ~ ag_df$clear_zoon_confirm, pch = 16, asp = 1,
       xlab = "Zoonotic outbreaks", ylab = "Spatial lag",
       main = paste("Moran's I -", name))
  M1 <- lm(inc.lag ~ ag_df$clear_zoon_confirm)
  abline(M1, col = "blue")
  dev.off()
  
  # Extract results
  morans_coef <- coef(M1)[2]
  robust_se <- sqrt(diag(vcovHC(M1, type = "HC4")))[2]
  mean_confirm <- mean(ag_df$clear_zoon_confirm, na.rm = TRUE)
  p_value <- coeftest(M1, vcov = vcovHC, type="HC4")[2,4]
  sig <- ifelse(p_value<0.05,"yes","No")
  
  results <- rbind(results, data.frame(
    subset_name = name,
    morans_I = morans_coef,
    morans_I_se = robust_se,
    mean_confirm = mean_confirm,
    p_value = p_value,
    significant=sig,
    stringsAsFactors = FALSE
  ))
}

# Show results
print(results)